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SUMMARY 


A method for analyzing shock coalescence which includes three-dimensional 
effects is developed. This method is based on an extension of the axisymmetric solu- 
tion. The asymmetric effects are introduced through an additional set of governing 
equations/ which are derived by taking the second circumferential derivative of the 
standard shock equations in the plane of symmetry. The coalescence method is consis- 
tent with and has been combined with a nonlinear sonic boom extrapolation program 
which is based on the method of characteristics. Though the two sets of governing 
shock equations are uncoupled, the flow equations, developed in the same manner, are 
weakly coupled through the first derivative of the circumferential velocity. Since 
the characteristic behind the shock is necessary to fix shock location, an iterative 
procedure between the axisymmetric and asymmetric equations is used. The extrapola- 
tion program, originally unable to handle shock coalescence, is now able to extrapo- 
late pressure signatures which include embedded shocks from an initial data line in 
the plane of symmetry at approximately one body length from the axis of the aircraft 
to the ground. 

Descriptions of the axisymmetric shock coalescence solution, the asymmetric 
shock coalescence solution, the method of incorporating these solutions into the 
extrapolation program, and the methods used to determine spatial derivatives needed 
in the coalescence solution are included. Results of the method are shown for a body 
of revolution at a small, positive angle of attack. The body was designed so that 
embedded shocks would be included in the near-field data. These results are corrpared 
with results from two widely used sonic boom propagation methods which are based on 
modified linear theory. Signatures are shown at several stages of the extrapolation 
process through the atmosphere in an attempt to understand the reasons for the dif- 
ferences in the results of the three methods. 
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U INTRODUCTION 


1.1 Purpose 

This paper presents an analysis and method of solution for shock coalescence 
with asymmetric effects. The analysis is consistent with and an extension of a sonic 
boom propagation theory developed at New York University (refs. 1 to 3) and has been 
incorporated into the coirputer program associated with that work. To provide infor- 
mation on the need for and the importance of this work, a brief review of sonic boom 
prediction, extrapolation, minimization, test methods, and some weaknesses of these 
methods is included in the next section. 


1.2 Review of Sonic Boom Work 

The sonic boom continues to be a major obstacle in the path of the development 
of an economically viable supersonic transport fleet. Though much progress has been 
made in understanding, predicting, and minimizing the sonic boom, the ultimate goal 
of producing an aircraft with sonic boom characteristics acceptable for overland 
flight remains elusive. 

Most currently used prediction methods are based on theories developed by Whit- 
ham (ref. 4) for a supersonic projectile, and by Walkden (ref. 5) , who extended the 
analysis to include lifting bodies. These theories combined with the supersonic area 
rule theory developed by Hayes (ref. 6) led to the generally accepted prediction 
methods described in some detail in reference 7. Originally, it was felt that all 
sonic boom pressure signatures would have attained the standard N-wave shape when 
they intersected the ground. Wbrk by McLean (ref. 8) indicated that for airplanes 
with extensive lifting surfaces this was not so. The "freezing" effects of the real 
atmosphere pointed out by Hayes (ref. 9) reinforced the idea that the signature may 
still retain midfield effects when it intersects the ground. A midfield signature 
retains effects of the airplane shape and thus offers this shaping as a possible 
avenue of minimizing the sonic boom. 

Using as a basis the sonic boom prediction method, the idea of a midfield signa- 
ture, and work by Jones (ref. 10) on the far-field minimum, Seebass and George devel- 
oped a minimization procedure (refs. 11 and 12) which predicted a minimizing equiva- 
lent area distribution for given flight conditions. Extensions to this procedure for 
the real atmosphere and nose bluntness relaxation are given in references 13 and 14. 
Working inversely, the designer can develop a model whose area matches the above 
boom-constrained distribution. This process was used in the design of three wing- 
body configurations, which were tested experimentally for sonic boom levels (refs. 15 
and 16) . Conclusions from that study indicated that the approach is valid, but that 
the boundary layer and wake effects must be considered in the design process. Two 
feasibility studies (refs. 17 and 18), in which systems, passenger loads, safety, and 
efficiency considerations were included, indicated that a low-boom aircraft was 
within the realm of possibility, but that extensive in-depth studies and trade-offs 
were needed. Certain characteristic features of low-boom aircraft seemed to evolve 
from the studies. These features included a large wing area, long root chord, posi- 
tive dihedral, twist, camber, and canards. 


As theoretical sonic boom studies have evolved, so have the accortpanying test 
methods. The earliest propagation methods were based entirely on far-field theory 
and thus required as input data a far-field N wave. In order for the pressure signa- 
ture to have traveled a sufficient number of body lengths to have attained its N wave 
form in the existing supersonic tunnels, model size was limited to approximately 
1 inch. With the advance in knowledge of the inportance of the midfield signature 
and the development of propagation methods which accepted this data, wind-tunnel 
model size became limited by the requirement that the data be axi symmetric. Since 
tests have shown that disturbances caused by volume are nearly axisymmetric at about 
three body lengths away and that those caused by the lifting forces have also lost 
much of their asymmetry, this lowered restraint generally allowed model size to 
increase to approximately 6 inches. Model builders can include many more features in 
a 6-inch model than in one which is only 1 inch. Even so, details such as twist, 
camber, and nacelles are still nearly inpossible to include on a 1/50-scale model 
with the accuracy needed. 

To increase the reliability of testing the latest low-boom transport designs, 
the models should be limited only by tunnel size and flow quality within the tunnel. 
This would allow model size to increase to roughly 2 feet in span and 3 feet in 
length. This size would allow one model to be used for both force and sonic boom 
tests. In addition to allowing improved accuracy in construction, this size also 
reduces the problems with the boundary layer and the sting support system (ref. 19). 

Propagation methods based on linearized theory become inaccurate for strong 
shocks where significant entropy changes occur and for high altitude and Mach numbers 
where the cumulative effects of second-order terms become significant (ref. 20) . The 
propagation method (refs. 1 and 2) developed at New York University (NYU) by Ferri, 
Siclari, and Ting eliminates many of these problems. Nonlinear terms have been 
retained, and altitude effects and entropy effects are included. The method of solu- 
tion is the modified method of characteristics (MMOC), in which step sizes on the 
order of several body lengths may be taken without destroying the accuracy of the 
program. The program includes nonlinear effects, which are important in the near 
field, but still does not account for nonaxi symmetric effects, which are also very 
prominent. A second program (ref. 3), developed by Ferri, Ting, and Lo, includes the 
asymmetric effects near the vertical line of symmetry - where sonic boom effects are 
strongest. This is acconplished by bringing in three-dimensional effects through 
derivatives in the circumferential direction which do not vanish in the vertical 
plane of symmetry. The method introduces the effects of asymmetry without the pro- 
hibitive difficulty of the full three-dimensional method of characteristics. 

As stated, the low-boom transport model will have a very complicated, rotational 
near field with embedded shocks from the wing-fuselage juncture, from engine pods, 
and from other parts of the aircraft. In the near field, these embedded shocks will 
not have coalesced into one shock. Though the NYU program accepts input at approxi- 
mately one body length, it cannot handle the coalescence of shocks, and thus, calcu- 
lations cease. 

Using the same approach as that used to include the three-dimensional effects in 
the MMOC program, this paper presents the analysis for axisymmetric and nonaxisymmet- 
ric shock coalescence. 
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1.3 Outline of Discussion 


As previously stated, the major contribution of this paper is the analysis and 
solution for shock coalescence with asymmetric effects, as described in section 3. 
However, because this method was developed to be used with the MMOC program, it is 
also necessary to describe in some detail the methods of that program, the method in 
which the shock coalescence is combined with the MMOC program, and the methods used 
to find additional information needed for the solution. 

To this end, a fairly detailed description of the equations and method of solu- 
tion of the MMOC program are given in section 2. Also included in that section is a 
description of the closure condition of the method. The limitations introduced into 
the outgoing characteristic to enable larger step sizes are also described. 

Section 3 includes the analysis and method of solution for shock coalescence - 
both the axisymmetric and asymmetric cases. Once the point of shock intersection is 
located, the axisymmetric solution is found and then the asymmetric solution. No 
iteration is needed, since the axisymmetric and asymmetric governing equations are 
not coupled. 

Spatial derivatives such as u^ and u r are necessary for the asymmetric 
coalescence solution, and these variables are not defined in the MMOC program. Sec- 
tion 4 describes the methods used to obtain their values. Used in section 4 is a 
method developed by Lin and Rubinov (ref. 21), which states that if the curvature of 
a shock is known, then the spatial derivatives behind the shock are also known. The 
system of equations for the derivatives becomes indeterminate for shocks of zero 
strength and, therefore, unsuitable for numerical computations with weak shocks. In 
the current analysis, it is assumed that the resultant shock is of finite strength, 
but the isentropic shock or expansion wave of the opposite family needed to equalize 
the system is so weak that continuity of the variables across the wave may be used. 
For shocks of zero strength, the curvature of the shock (or that of the centerline of 
an expansion fan) differs from that of the characteristic ahead of and behind the 
shock because of the jump of the curvature of the streamline, which is continuous in 
slope. To generalize the system for shocks of any strength, it is necessary to 
regroup the equations in such a way that the common factor which vanishes for zero 
shock strength cancels analytically. 

The procedure for combining the shock coalescence with the MMOC program is 
described in section 5. This section includes a description of how the intersection 
point is found, how data immediately in front of the intersection point are deter- 
mined, the solution, and then how the solution is recombined with the characteristic 
network. 

Section 6 describes results which were obtained for a lifting body of revolution 
with embedded shock coalescence during propagation. It is beyond the scope of this 
paper to conduct a parametric study with many different bodies. Rather, the emphasis 
is placed on the development and incorporation of the method of shock coalescence 
with the MMOC program, and a simple body of revolution is used to illustrate the 
results. Extrapolated data from the MMOC program are compared with results from two 
well-known sonic boom prediction methods which are based on modified linear theory. 

This work was done to satisfy in part the requirements for the degree of Doctor 
of Science at George Washington University, February 1983. 
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2. THEORETICAL BACKGROUND FOR MMOC SONIC BOOM PROPAGATION METHOD 


2*1 Governing Equations and Closure Conditions 

The analysis developed in this paper is an extension of the methods developed in 
references 1 to 3. Governing equations for the steady-state, inviscid flow field 
including gravitational terms and variation in total energy are (ref. 22) 

Momentum 


V 


VV = - + G 

P 


( 2.1 ) 


Continuity 

pV • V + $ • Vp = 0 


( 2 . 2 ) 


Energy 

V • VH = 0 (2.3) 

State 

p = pRT (2.4) 


where the total enthalpy H is defined by 


2 

H = RT + 3— + gz (2.5) 

y - 1 2 

In these equations, g is the acceleration due to gravity, G is the body force 
vector per unit mass (g = -g Vz, where z is the vertical coordinate), R is the 
gas constant, and q is the fluid speed (q = |v|). 

It is convenient for computational purposes to transform the set of governing 
equations. The Gibbs relation (ref. 23) is introduced in the form 

T VS = - VT - — Vp (2.6) 

y - 1 p 

where^ S is the fluid entropy per unit mass. Use of equation (2.6) and the relation 
V. W=v(q/2) + (V x v) x V allows pressure and density to be eliminated from the 
momentum and continuity equations, which become 

(v X v) x V = T VS - VH (2.7) 


4 



and 


V 



$ • S 


2 

a 


( 2 . 8 ) 


2 

in which the speed of sound a is given by a * = yRT. It is noted here for later 
use that in view of equation (2.3), the component of equation (2.7) along a stream- 
line is 


V • VS = 0 


(2.9) 


Ahead of the bow shock, the flow is assumed to be uniform with velocity U • 

CO 

The temperature T (z) and the pressure p (z) are given functions of the altitude 
z which satisfy the condition of hydrostatic equilibrium 

T (z) dp (z) 

oo CO g 

p (z ) dz ~ ~ R 

00 

Jump conditions for the system consisting of equations (2.3) to (2.8) are discussed 
in detail later. However, it is pointed out here that the total enthalpy remains 
constant across a stationary shock wave (ref. 23). This, together with equa- 
tion (2.3) implies that H is invariant along streamlines 

2 u ® 

H = r RT + + gz = - RT f z 1 + -r— + gz (2.10) 

Y ~ 1 2 13 Y - 1 oo v cdj 2 00 

where the notation z is introduced to denote the elevation of a particular stream- 

00 

line in the uniform flow ahead of the bow shock. In the numerical computations 
described later, the function z^(x r y,z) was actually calculated instead of H 
itself, by replacing equation (2.3) with V • Vz =0. Once z or H and the 
velocity are known, then T can be considered as known from equation (2.10). Thus, 
the system consisting of equations (2.3), (2.7), and (2.8) can be regarded as a sys- 

tem of five equations for five unknowns: u, v, w, S, and H. 

For an airplane flying at constant altitude Y, it is convenient to introduce 
cylindrical coordinates x, r, <|>, shown in sketch A, with velocity components u, 
v, w where 


z = Y - r cos (\) 


and 


A 



5 



T 

Y 

1 


T 

d 

z 

1 



Sketch A 

In this coordinate system, the five fundamental equations (2.3), (2.7), and (2.8) are 

expressed as 


dH dH w dH 

u — + v — + — — 

dx or r o<\> 


0 


( 2.1 1 ) 


du 5v ;v J_ dw 
dx + dr + r + r d4> 



dx V dr 


w 

r 


d \( u 2 + v 2 + w 2 

a <!7 \ 2 


+ gz 


( 2 . 12 ) 


T 


as 

5x 



t as 

r dc); 


9H _ /du dv\ h du_ dw\ 
dx ~ v \pr ” 6xJ + w \r dc|j " dxj 


dH _ f dv 8u\ /l_ dv dw _ w\ 

5r ” u \px ” dr/ + yr d(p “ dr r/ 


J_ J_ du\ / dw J_ dv 

r 5([) ” U \dx ~ r 5ci»y + ydr r 5cp 


(2.13) 


(2.14) 


(2.15) 


In the vertical plane of symmetry, the circumferential velocity w vanishes as do 
the first derivatives of all variables except w with respect to the circumferential 
direction. Thus in the plane (\> = 0, equations (2.11) to (2.14) become 




0 


(2.16) 


J\ au 

a 2 ) 5x 


v | 5v uv/du dv \ v _1_ d_w cjv _ 

a 2 J dr a 2ydr + dx/ + r + r d(|; + ^2 " 


dS dH 
dx “ dx 



(2.17) 


(2.18) 
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(2.19) 


as 5H / Bv _ Bu 

dr ~ &r u \Bx ar 


Equation (2.15) is identically satisfied in the plane (\> - 0. 
Equations (2.18) and (2.19) are combined to give 

Bv au 1 / as bh \ 

ax Br “ q \ Bn Bn I 


( 2 . 20 ) 


where q B/Bn = u(B/Br) - v(B/Bx), and B/Bn is the derivative normal to the stream- 
line. Equation (2.20) represents one scalar component of the original vector equa- 
tion (2.7). The second independent relation will be taken to be equation (2.9), 
which in the plane of symmetry becomes 


as as 

u Bx + v Br ~ 


0 


( 2 . 21 ) 


Equations (2.16), (2.17), (2.20), and (2.21) now form a set of four equations 
with five unknowns: u, v, S, H, and Bw/B^. For the axi symmetric approximation, 

it is assumed that the term (1/r)(Bw/3<J ;) is negligible at large values of r. Thus, 
a closed system of equations is attained. If this assumption is not made, then 
another method for closing the system must be found. 

Since equation (2.15) is an identity in the plane 4> = 0, its first derivative 
with respect to c|> is taken. This process introduces an additional equation, but it 
also introduces several new unknowns. The second derivative with respect to <\) is 
taken of equations (2.11) to (2.14) and (2.9). The derivatives of equations (2.13) 
and (2.14) are combined by multiplying (2.13) by u/q and subtracting the product of 
v/q and equation (2.14). This process introduces five new equations and five new 
unknowns: B^u/Bc [ , B^v/Bcp^, B^H/Bc^, B^S/B(|>^, and B^w/B4>^. Ihe entire system 

now consists of 9 equations and 10 unknowns. 

The closure assumption introduced in references 1 to 3 is 


B 3 w Bw 

a, 3 * ' a * 


( 2 . 22 ) 


This assumption is justified on the basis that the primary contribution to the cross 
flow can be described by a function of the form 


w - w(x,r) sin <|> 


which is a valid approximation describing asymmetries due to lift or small asymme- 
tries which are present at approximately one to two body lengths from an aircraft. 
With the closure condition, the system reduces to nine equations and nine unknowns. 
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The last five equations formed by taking derivatives in the circumferential 
direction are listed below where the closure condition, equation (2,22), has already 
been applied. For convenience, the following notation is introduced: 


u(x ' r) l 1= 0 


U = 


Uc 


5u 


d 2 u 


dw 


:0 - 8X 

6=0 

„ 5^ 2 

cp=0 


6=0 

X U a 

Uoo 



Equations (2.11) to (2.15) yield 


uH + uH + vH +vli + 2w h = 0 
x x r r r 


(2.23) 


u + v — 

X r r2 


"2 2 

uu + UVfU + V ] + V V 

x \ r xj r 


Z + 2t + 

r r a 2 


u ^2uu x + vu 


v ^2vv^ 


+ UV X + uu 


r + vv x) 

t) + g(v - v) 


+ 2w(vw^ + uw x + g^ + ~~ fw + uu 


V r 


V 


A 2^ ~ 2 (- v - w\1 

+ w / " a V u x+7+v r +-J 


(2.24) 



TS + H + uu + vv = urw x + vrw r + vw 


uS + uS + vS + vS + 2w s = 0 
x x r r r 


(2.25) 


(2.26) 


(2.27) 


Equation (2.25) results from the combination of the derivatives of equations 
and (2.14). The operator d/dn is defined beneath equation (2.20). Bjua- 


(2.13) 
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tions (2.16), (2.17), (2.20), (2.21), and (2.23) to ( 2 .J2 7) orm a system of nine 
equations and nine unknowns: u, v, w, H, S, u, v, H, and S. The first set, 
equations (2.16), (2.17), (2.20), and (2.21), are coupled weakly to the second set, 

equations (2.23) to (2.27), through the term ( 1/r) (dw/d(J>) in equation (2.17). 

The method of characteristics (ref. 24) is applied to equations (2.17) 
and (2.20) and results in the following characteristic equations: 


cot 


ae 
p q 


52^(5 ds 



sin u dr 
sin(0 ± \x) 


sin 0 




0 


(2.28) 


These equations are valid along the characteristic directions 


= tan( 0 ± ji) (2.29) 

dx v 

| -2 -2 - - 1 - - 

where q =y u + v , 0 (the direction of the flow) = tan (v/u) , and 

\i (the Mach angle) = sin ^(l/M). 

A second set of conpatibility equations are derived from equations (2.24) and 
(2.25). Equation (2.24) is expressed as 



(2.30) 


where 9* stands for the right-hand side of equation (2.2 4) rewritten in the form 



2w /- 


(ts - 


H + 2uu + 2vv - 


— A 

vw 


-2\ 
" ) 


- ~ * - a 

+ g I v - v + 2w + v — 


- »~ 2 » 2 b 
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in which 5/5 a is the derivative along the streamline defined by 

q5/5a = u(S/5x) + v(5/5r). Then equations (2.25) and (2.30) yield a second set of 
coirpatibility equations 


du + tan (3 + \i) dv « sxn_jx — r (±J cos \i - gr sin \i) 

sin( 0 ± \i) cos ( 0 + p.) 

along the same characteristic directions 


(2.31) 


= tan(0 • \i) (2.32) 

dx 

where J represents the right-hand side of equation (2.25). 

Equations (2.28), (2.16), and (2.21) provide solutions for u, v, S, and H 

for a given value of w, and equations (2.31), (2.23), (2.26), and (2.27) yield 

u, v, S, H, and w. Solutions are obtained by simultaneously solving for the 
velocity components at the characteristic intersections and for enthalpy and entropy 
along streamlines. 


2.2 Shock Conditions 

The treatment of the bow shock and all embedded shocks is approached in the same 
manner. When the complete set of five shock equations are limited to the plane 
(\) = 0, the (^-momentum equation vanishes indentically . There are four remaining equa- 
tions to solve for the shock shape f^ and u, v, p, and T behind the shock. 

The compatibility condition along the C + characteristic behind the shock provides 
the additional equation needed to close the system for the axisymmetric case. The 
equations needed to determine u, v, w, T, and p are derived by taking the first 
derivative with respect to <\> of the ^-momentum equation across the shock and the 
second derivative with respect to <\> of the remaining shock equations. (The proce- 
dure for deriving these shock equations is described in detail in Section 3.) This 
process introduces six unknowns, u, v, w, p, T, and f , and five new equations. 
The compatibility equation for the second derivative quantities, equation (2.31), is 
used to close the system behind the shock. 


2.3 Method of Solution 

A 

Because the two sets of equations are coupled only weakly through w, the method 
of solution is iterative. A value of w is assumed^ and the system_consis ting of 
equations (2.16), (2.21), and (2.28) is solved for u, v, H,^and^ S.^ This 
information is used to solve the second set of equations for u, v, H, S, and w. 
The entire procedure is repeated until the value of w converges. 

The program is an extrapolation program, and it requires initial input approxi- 
mately one body length from the axis of the aircraft. These data must be provided 
either by experimental methods or by programs which calculate the entire flow field 
about an aircraft. From the first data points on the initial line, a new bow shock 
point is located. Calculations proceed from this point along the C~ characteristic 
until the initial data line is intersected. Thus, a new data line is created. This 


1 0 



process is repeated until finally the new data line is a C“ characteristic which 
terminates at the last data point on the initial line. The next C“ characteristic 
is then started, and calculations continue to ground intersection. (See sketch B.) 


t 



2.4 Far-Field Modifications 

Because of the extreme propagation distances necessary with this method and thus 
the possibility of accumulated error in many steps taken, the outgoing 0 + \x char- 
acteristics are modified according to Whitham theory to allow for a gradual increase 
of step size in C + to the order of several body lengths as r increases. Within 
this modification, the characteristics are no longer considered to be straight but 
are curved in a manner consistent with Whitham theory. The modification is applied 
at each step in the 0 + \i characteristic and along the shock waves. 

Whitham (ref. 4) gives the following relationship for the slope of the charac- 
teristic 

k F 

dx _ oo 

dr ~ Pen 1/2 

2r 


where 


p. 



k_ = 


_ y + 1 


M 


F 


r 


1/2 Vp 
P 


The exact characteristic is modified by 


dx 

dr 


k F 


2r 


1/2 


- cot( 0 + |i) - |3 


k F 
00 


2r 


1/2 


The characteristic equation is then rearranged so that the left side is an exact 
differential of the invariant under Whitham theory, and the right side is then an 
order of magnitude smaller than that in a standard characteristic equation 



( 2 . 33 ) 


d ^ x " Pco r + k oo Fr1/2 ^ = jc° t{0 + + 

Thus, 

x = 6 r - k Fr 1//2 + /e dr + A 

r CO oo — c 


where 

k F 

E = cot(0 + \l) - + — jj 2 

and A c is a constant of integration. E is smaller than cot(0 + |i) and thus 
allows the step size to gradually increase to several body lengths in the far field. 
A similar procedure is used for the bow shock and for embedded shocks. A more 
detailed description of this procedure is found in reference 1. It is determined in 
reference 1 that errors of less than 1 percent are present in ground signatures even 
when the radial step sizes are allowed to increase gradually to several body lengths 
in the far field. 



3. COALESCENCE OF SHOCKS 
3.1 Axisymmetric Shock Coalescence 

3.1.1 Governing equations .- Once removed from the axis of the body, the treat- 
ment of the shock interaction for the axisymmetric case becomes the same as for the 
two-dimensional case. Thus, when the shock angle is defined, the oblique shock equa- 
tions may be used. 

Incoming shocks f^(r) and f 2 (r) are gi ver h and it is assumed that the point 
of intersection p Q has already been determined. (See sketch C.) Theory tells us 
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that from a shock-shock interaction of this nature, we may expect a resultant shock 
f 3 (r) of the same family, and a slipstream h created because of the difference in 
entropy of the flow passing through one shock and that passing through two shocks. 
Because the pressure and flow direction must be equal at the slipstream, a weak 
(isentropic) expansion or shock f 4 (r) of the opposite family may be needed to allow 
both conditions to be met. 


In this situation, the unknowns are u 4 , v^, p 4 , T^, v^, u^, p^, T^, (3^, 

and p 4 , where p is the shock angle with respect to the flow ahead. The standard 
oblique shock equations are 


JLtH 

J 5,n 



1/n 

+ 5 
n 



(3.1) 


(3.2) 


v 


1 


r t 


V 


5, t 


(3.3) 


T 5 K.n - 0( M 1.n + 5 ) 


36M 


1 ,n 


The conditions holding across the slipstream h(x) are 


(3.4) 


and 



P„ = P. 


(3.5) 

(3.6) 


and across the isentropic wave in regions 3 and 4, we require 


T 


4,T 


T 


3,T 


(3.7) 


P 4,T P 3, T 


(3.8) 



(3.9) 


(3.10) 
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The system consisting of equations (3,1) to (3,10) is closed with 10 equations and 
10 unknowns. All values in regions 1 and 3 are known, 

3.1,2 Method of solution ,- The solution process is begun by assuming an initial 
value of (3^; the guess used is 


h - <p, + iy /2 

Using equations (3,1) to (3.4), u^, p^, v^, and T 5 are calculated. From equa- 
tions (3.5) and (3.6), values for 0 and p are obtained. Because total pressure 

4 4 

and total temperature in regions 3 and 4 are equal and because all information in 
region 3 is known, equation (3.9) may be applied in both regions to calculate 



and equation (3.10) can be used to calculate 


T 3 |i.O ♦ (m 2 3 /5)_ 

1 .0 + I M^/5 ) 

V 4/ / 


Since this is a perfect gas 


= ^yRT 


U 4 = 


a.M . 
4 4 


yl + tan' 


and 


v = u tan 0 
4 4 4 


(3.12) 


(3.13) 

(3.14) 


(3.15) 


The Prandtl-Meyer angle in both regions is 


v 


4 




- tan 




(3.16) 


v 


3 





tan 


-1 




(3.17) 
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The angle of the weak wave f 4 is calculated using the standard pressure ratio 
equation 



It has already been determined that the wave f 4 must be of the opposite fam- 
ily. If it is an expansion, then p^ < p 3 , and therefore, > M 3 and v 4 > 

Also, because it is of the opposite family, then © 4 > 0^. Similarly, if it is a 
shock, then p 4 > p 3 , M 4 < M^, < v^, and 0 4 < 0^. The error in the original 

guess is assessed by calculating 

Error(N) = (v. - vj + (0 O - 0_) (3.19) 

v 4 3' ^ 3 4 y 

To iterate, a new value of (3 is assumed, and all calculations are repeated. After 
the second iteration, a straight line extrapolation is used to guess a new value of 
(3 3 ( N) 


P 3 (n) 


Error ( N- 2 ) p ( N- 1 ) - Error(N-l) p (N-2) 
Error(N-2) - Error(N-l) 


(3.20) 


where N represents successive iterations after the second. Iterations are 
continued until 


Error(N) < e 


where e is some predetermined error criterion. Generally, convergence occurs 
within five or six iterations even when the error criterion e is on the order of 
1CT 10 for weak shocks. Once the angles P 3 and p 4 have been determined, the 
entropies S 4 and S 5 may be calculated using standard shock equations. The slopes 
of the shocks and slipstream are 


3 r tan( p 3 + 0^ ) 

f = 

4 r tan <P 4 + e 3 ) 

and 



(3.21) 

(3.22) 

(3.23) 
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3.2 Asymmetric Shock Coalescence 


3.2.1 Derivation of shock equations. - The shock conditions across the resultant 
embedded shock x = f^lr,^) are now developed. The region ahead of the shock is 
denoted by the subscript 1 and that immediately behind the shock by the subscript 5. 
The conservation conditions across the shock are 


Mass 


Pi 



V 1 f 3 



= m 




(3.24) 


Momentum 


mu + P 1 = mu 5 + p 5 


V 3 + v , ’ V 3 + v 5 

r r 


f 3 . 

(j> (j) 

U, + W. = U c ^ + w_ 

1 r 1 5 r 5 


(3.25) 

(3.26) 


(3.27) 


Energy 


2 2 2 2 1* T , 

u„ + v. + w. + r 

1 1 1 y - 1 


2 2 2 2 ^ RT 5 

U 5 + V 5 + W 5 + 7~^T 


(3.28) 


To avoid difficulty as the shock strength becomes weak, these equations were 
rearranged in reference 3 to yield the following three governing equations: 


( u 5 - u i )f 3 + v 5 - Vl = 0 


f 3 

(u 5 - ) “7^ + w 5 " w 1 =0 


and 


(y + 1 ) Q-| Q 5 - (Y - 1)Q 1 = 


2 Y RT l\1 + f 3 + ~2 Li 

r r 


(3.29) 


(3.30) 


(3.31) 
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where 


w. f 
x 3 


Q. = u . - v.f ^ 
*1 1 l 3 


<1* 


These three equations are used to provide solutions for u, v, and f r in the plane 
of symmetry. To determine p and T, equations (3.2 5) and (3.28) must be used. 

The asymmetric shock conditions are developed when each of these equations is 
differentiated twice along the shock surface with respect to cp, and then (p is set 
equal to 0 to reflect conditions in the plane of symmetry. In the plane <|> = 0, the 
first derivatives with respect to cp of all variables and shock surfaces vanish with 
the exception of w.. A detailed description of this derivation will be given for 
shock equation (3.29). 

Let 


A - (“5 - u , > f 3 + V 5 ‘ V 1 * 0 


The derivative along the shock surface x 

5f 3 d_ + a 


f^(r,tp) with respect to cp is 


dcp d(p 5x dcp 


The first derivative with respect to cp is 


_5A_5A _5A 

Q cJj — 3x f 3 + 5^ A x f 3 + A <|> 

The second derivative then yields 


d 2 A 5 


5<|; 


(A f , + A. I f- + -r— /A f , + A . \ 
\ / <\> \ / 


Note immediately that the first term in this expression will vanish in the plane 
cp = 0, since the factor f^ vanishes. Thus, this term will not be expanded 

further. Working only with the second term, then 


* (Vs„ * \) 


= A f + A f + A 

^ 3 <(, x 3 4>4> W 


Again, the first term of this expression also vanishes in the plane cp 
Therefore, in the plane cp = 0 


= 0. 


-Mr = A f + A 

a 1 2 x 3 I I 44 
dcjj (Jjcp 
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Expanding A^ and A^ yields 


A =» u f + u f - u f - u f + v - v 
x53 53 13 13 5 1 

x r rx x r rx x x 


A = u f + u_f - u f - u f + v - v 
(\) 5.3 53 , 1,3 13 5 1 

^ r r<\> r r<\> (|> <|> 


A, , = u r f + u f + u r f + u_f - u f 

(btb 5 . . 3 5,3, 5,3. 53,. 1 . , 3 

^ 4>cp r cp r4> cp rrcp r*(p cp cpcp r 


+ u f - u f - u f + V - V 

1 cp 3 rc|; cp rep r4><J> cpcp cpcp 


In the plane cp = 0 


A 

x 


u 

1 



+ V 


V 

1 


X 


and 


<14 


= u 


44 


U 5 f 3 


- U f - U f + V - V 

rep cp cpcp r r cp cp cpcp cpcp 


Therefore 


5 2 A | 

54> 2 


4^=0 


K f 3 ’“I f 3 +V 5 - V 1 ) £ 3,, + U S £ 3 

\ x r x r x x/ cpcp c|4 r 


+ U 5 f 3 ,, - V , f 3 - U 1 f 3 .. 

cpcp r r4>4> 


The bar (” ) and tilde (~) over symbols denote values in the vertical plane cp == 0 • 
Thus, the asymmetric version of equation (3,29) becomes 


("5 x ? 3 r - 


^x ? 3 r + ^5 X - v lx )?3 


(V 3r 


u 5 f 3 r " U 1 ^ 3 r “ u 1 f 3 r + v 5 " V 1 


)■ 


(3.32) 


By a similar method, equation (3.30) becomes 

£ 3(^1 - u 5 ) 
w r = + w. 


(3.33) 
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ll 


where 


dw 

1 

* d<i> 

f 

0 


and equation (3.31) becomes 


27R 


(V3 + ^ + f 3r) + 4 ^ RT ll f , f , + 


f 3 


3 r 3 r r 2 


(<2iQ 5 + Ql25 ij - 


= <Y + 1) Q 1 Q 5 + Q 1 Q 5 - 2(y - DQ^! 


(3.34) 


where 


rv* — ~ — — « fv •_ _ «j 2 W4 

Q. = u f - V. f f + u. - V.f - v. f - — yT- 
i x ^ i x 3 j- 3 i i 3 j~ i 3 ^* 


and 


2i 


u. - V. f 

1 1 

r 


(i = 1 f 5) 


3.2*2 Asymmetric analysis *- The three-dimensional geometry of the shock coales- 
cence is described by sketch D, and the geometry in the plane <|> = 0, by sketch E. 


Resultant shock 





Sketch D 


19 



Sketch E 

Incoming shock waves x = f ^ (r ^ cp ) and x = coalesce to form the following 

system: resultant shock wave x = f (r, 4 >), slipstream r = h(x,(|;), and a weak 

(isentropic) wave of the opposite family x = f (r , cp) • This results in five surfaces 
intersecting along the line x = X(<p) # r ~ R ( • Taking the second derivative along 
each of the intersecting surfaces with respect to c|> and setting <\> equal to 0 
yields 


X = f j R + f j (j = 1, 2,3,4) (3.35) 

R = h x X + h (3.36) 


Along the slipstream surface r = h(x,(Jj), pressures and flow direction on either 
side are the same. Thus, 



-> a 

V 4 • n = 0 


-> * 

V 5 • n = 0 


where n is the normal to the slipstream defined by 


n = j -hi 
J x 



The equations for flow direction at the surface h become 


-u. h + v. 
IX 1 


w, 


-4-0 


1 r 


(i = 4, 5) 
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Second derivatives along the slipstream surface with respect to ()> are calculated by 
using 



5_ 
44 dr 


+ 



and thus the pressure equation yields 


P4 r ^ + P 4 = P5 r ^ + P5 


(3.37) 


and the flow equations 


a ~ 

2 wd h 


-u hh + v.h-u.h - u.h + v, - 
a r 


4 x J. J. X I V . VX jli. w. I V J 

r x 4 r 4x 4x 4 r 


= 0 


(3.38a) 


__ _ ~ ~ n _ m rw 2 w^h 

“U _ h h + v_ h - u_h - u^h + v_ - : — = 0 

5 r x 5 r 5x 5x 5 r 


(3.38b) 


The three shock equations for the asymmetric quantities (eqs. (3.32) to (3.34)) were 
derived in the previous section and are valid from regions 1 to 5 across shock f-. 
To calculate T behind the shock, the original energy shock equation (eq. (3.28)) 
may be used. Taking the second derivative of this equation with respect to and 

setting c|> equal to 0 yields 


2 u,- u 


5, ' 2 U 1"1 X + - 2 v 1 v 1 x * Y 


(% - 


To 


I — A# — /V 

1 2 u u — 2 u u 
L 5 5 11 

(f^)(*5 ' *,)] ’ ° 


^ 2 ^ 2 — ~ - « 
+ 2w 5 - 2w i + 2v 5 v 5 “ 2v 1 v 1 


compute S 5 across the shock f^, the method in reference 2 is used 


S 5 " S 1 MS - S 1 

X 


(®5 ” ®1 V^3 + 2 

\ x x/ dM. 

' 1 ,n 


AS sh 2 . 2 „ 

— - M 1 sin p 


(3.39) 


(3.40) 
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where AS 


is the jump in entropy across the shock 
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one gets 



This equation can then be used to get 


p in region 5, or 



(3.41) 


For the examples shown in this paper, the surface x = f^(r,c};) is an extremely weak 
(isentropic) shock or expansion wave across which the variables are nearly continu- 
ous. We therefore introduce the approximations that 


u 3 " u 4 
v 3 = v 4 

P 3 = P 4 
s 3 = s 4 
T 3 = T 4 

A A 

W 3 = W 4 

Taking the second derivatives of the above relations along the weak surface f^ 
yields the following five equations: 


^3j/4 + u 3 = u 4x f 4 + u 4 

v„ £ + v = v „ f + v 
3 X 4 3 4 X 4 4 

P3j/ 4 + P3 = P4 X ^4 + P4 

s 3 x f 4 + S 3 = S 4 / 4 + S 4 


T 3 x f4 + T 3 


T4/4 + t 4 


(3.42) 

(3.43) 

(3.44) 

(3.45) 

(3.46) 
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and the following equation is retained for w: 


W 3 = w 4 (3.47) 

The unknowns in^the asymmetric^system are X, ^ R, f 3 , f 4 , h, f 3 , h x , 

u 4 , v 4 , w 4 , T 4 , S 4 , p 4 , U 5 , V 5 , w 5 , T 5 , S 5 , and p 5 . The equations valid 
in the system are the five intersection equations (eqs. (3.35) and (3.36)), the two 
slipstream equations (eqs. (3.37) and (3.38a)), the six shock equations (eqs. (3.32) 
to (3.34) and (3.39) to (3.41)), and the six continuity equations from regions 3 to 4 
(eqs. (3.42) to (3.47)). There are 19 unknowns and 19 equations for the system. The 
system is closed. The slipstream flow equation (eq. (3.38b)) has been discarded 
because of the continuity approximations made across f 4 » It is found that the 
approximation of continuity described above yields a rapidly converging solution 
whenever | (p 4 /P 3 ) - 1 | < 10 “^. 

When | C p 4 / p 3 ) - l| > 10 ~ 5 , it: is necessary to treat f 4 as a finite shock or 
expansion. If f 4 is a shock, then the complete set of shock equations (eqs. (3.32) 
to (3.34) and (3.39) to (3.41)) instead of equations (3.42) to (3.47), are applied 
from regions 3 to 4, equation (3.38b) is used, and a new unknown f 4 ^ is introduced. 

The system then contains 20 equations and 20 unknowns. If f 4 is a finite expansion 
wave, the appropriate isentropic expansion equations are used. 

3.2.3 Asymmetric solution process. - For the incoming shocks f^ and t 2 ' t5ie 
values of fl r ' ^2 r ' f 1 r ^2' fi r # and ^2 r are a l rea< 3y known from the MMOC 
program. Equation (3.35) can be applied to shock surfaces 1 and 2 to solve for 


R 


h - 

fl r ^2 r 


(3.48) 


and 

X = f 1r R + f 1 (3.49) 

Knowing X and R, equations (3.35) and (3.36) may be applied to shock surfaces 3 
and 4 to solve for 


fi = X - f i R (j = 3, 4) 

J Jr 


(3.50) 


and 


h = R - h x X 


(3.51) 
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where f , f , and h are known from the axisymmetric solution. 

3 4 x 

r r 

thing is known in region 1 , all bar (") quantities are assumed known, 
known, then equation (3.3 3) may be immediately used to solve for W 5 ; 


Since every 

and ?3 is now 
i.e. , 


w_ = 


g 3 ( U 1 ~ "5) 


+ W^ 


(3.52) 


Similarly, since all quantities in region 3 are known, equations (3.42) to (3.46) may 
be solved for the tilde ( ~ ) variables in region 4. (If equation (3.47) is tested and 
the shock equations are needed at f Ar a double iteration procedure is used to deter- 
mine an additional unknown f 4 .) These quantities are 


u 4 = U 3 X ^4 + u 3 “ U 4 X ^4 (3.53) 

^4 = v 3x f 4 + v 3 - v 4x f 4 (3.54) 

P4 = P3 x f ~4 + P3 " P4 X ^4 (3.55) 

T ~ 4 = T 3 x f 4 + T 3 - T 4x f 4 (3.56) 

s 4 = S 3x f 4 + s 3 - S 4x ? 4 (3.57) 


To solve for the other tilde (~) variables in region 5, an initial guess is made 
for the value of £3 


f 3 = V2 f, + f 2 ' 

r V r r j 


Equations (3.3 2) and (3.34) are then used to solve for u^ and v^ 
^ U 1 - l 2 + l 3 y)^3 r + L 1 


v 5 = 


1 H- f§ 
J r 


(3.58) 


where 





L 1 - 



L 

3 


= v f 
5 3 


r 



- v f 


5 

x 


3 

r 



(3.59) 
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L 


2 


B 


5 



Li 


v i 


+ 





+ Vi 
'x 



*3 


B 


5 


B 4 ' B 3 + B ! + B 2 


B 4 = 2(y - DQtQ! 


B3 


B 


2 


(y + i )Q iQ5 



B-i = 2y 1 + 



Q 5 = U 5 " v 5 f 3 


2 1 = U 1 - V 1 f 3 


Q 1 = U 1 f 3 ‘ V 1 f 3 f 3 + U 1 " V 1 f 3 - V 1 f 3 


2 *1 ? 3 


x x r 


r r 

Using equation (3.39), T 5 may be obtained as 


T = T - 
5 1 


(VV * *,)] 


(3.60) 


where 


4. = 


U 5 U 5 


Vi 


V 5 V 5 


V 1 V 1 


+ V)(\ ' \) 


0 


2 


u u 
5 5 


Vi 


+ w 


-2 

w 

1 


V 5 V 5 


v v 
1 1 
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Equation (3.40) is used to solve for S 5 , and equation (3.41), for p 5 . The error in 
the iteration is evaluated by using equation (3.37) 


Error(N) 




P 


5 


A new guess for 



is made, for example. 


f 3r (2) = 1.1f 3r (1) 

and all variables from V 5 forward are recalculated with the new f 3 . From the 
values ? 3 r ( 1 ), error ( 1 ), £ 3 ^( 2 ), and error ( 2 ), a straight line approximation is 

used to determine the new fo 

J r 

f 3 r (N-2) Error (N-1 ) - f 3 r (N- 1 ) Error(N-2) 

^3 r ^ N ^ Error (N-1) - Error(N-2) 


The process is repeated until 


Error ( N) < e 


where e is a predetermined error criterion. For the examples shown, even if e is 
1CT 10 , convergence occurs within five iterations. After convergence, equa- 
tion (3.38b) is used to calculate h x * 


h = — -u c h h + v c h - Ur-li + v c - 


2w hi 
b 


5 X 
r 


5 11 

x 


(3.61) 


If the shock equations are used at f^, then equation (3.38a) is used to check 
convergence 


-u 4 h x h + v 4 h - u 4 h x u 4 h x + v 4 - 


2w h 
4 


0 


The solution process for the asymmetric equations is summarized in the flow 
chart shown on the next page. The solid lines on this chart represent the case in 
which the wave f 4 is so weak that continuity may be assumed for u, v, p, T, S, 
and w. Hie dashed lines are used to represent the double iteration procedure which 
must be used when the shock equations are used at f 4 . In that case. 
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Use eqs. (3 

.48 

to 3.51) 

to calculate X, 

R, 

f 3 , fj, ana h 


Use eq. (3.52) 
to solve for w r 


Use eq. (3.52) 
modified to 
solve for w. 


Use 

eqs . 

(3.53 to 3.57) 


to 

calculate 

V 

v 4 ' 

P 4 , T 4 , and S 4 


Use shock eqs. (3.58 to 3.60), 
(3.40), and (3.41) modified to 
calculate u 4 , p 4 , ? 4 , and £ 

between regions 3 and 4 


Guess f. 


Use eqs. (3.58 to 
3.60), (3.40), and (3.41) 

to evaluate 


V V 


T 5 , S 5 , and p 5 


Guess 
new f _ 


eq. (3.37) 


Calculate h 

x 

with eq. (3. 61] 


Solution 

complete 


Solution 

complete 


4 

weak? 


eq. ( 3 . 38a) 
\ = 0 ? ^ 














equations (3,53) to (3.57) must be replaced by equations comparable to equations 
(3.52) , (3.58) to ( 3 .60 ) , (3.40), and (3.41)) to solve from regions 3 to 4. 


4. SPATIAL DERIVATIVES NEEDED FOR ASYMMETRIC SOLUTION 
4.1 Spatial Derivatives in Regions 4 and 5 

The system of equations for nonaxi symmetric shock coalescence (eqs* (3.24) 
to (3.47)) is a closed system with the assumption that all spatial derivatives of 
u, v, p, and T are known. The method used to obtain these derivatives in 
regions 4 and 5 is now explained. 

In region 5, the following four governing flow equations are valid: 


x- Momentum 


p U U * p V U + TP - 0 
x r x 


(4.1) 


r-Momentum 


Ws + 5 sVs * ^ 5 5 - 

x r r 


(4.2) 


Continuity 


P 5 T 5 U 5 + P 5 T 5 V 5 ‘ U 5 P 5 T 5 " V 5 P 5 T 5 


+ U 5 T 5 P 5 + V 5 T 5 P 5 


P 5 T 5 


5 + h) 


(4.3) 


Energy 


+ u_, v v + u V u c + v^v c + — u T + v_T. = gv_ 

555 555 555 55 y-1 55 v - 1 55 5 

x x r r 1 x 1 r 


(4.4) 


Because of the weakness of the wave between regions 3 and 4 and the continuity 
assumption made, ^4 rr 3oes not appear as an unknown in the system. We therefore 

apply only three flow equations in region 4 - the two momentum equations and the 


29 



continuity equation. The energy equation is used only as a check for the continuity 
assumption. The equations in region 4 are 


P-U.U + Pa V a U a + T aPa = 0 

4 4 4 4 4 4 4 4 

x r x 


P « u * v . + p.v.v. + RT p = -gp„ 

4 4 4 4 4 4 4 4 4 


P 4 T 4 U 4 + P 4 T 4 V 4 - U 4 P 4 T 4 ‘ V 4 P 4 T 4 * U 4 T 4 P 4 
x r x r 3 


- - - P 4 T 4 (* - \ 

+ V 4 T 4 P 4 = ' r \ 4 + V 4 ) 
r v 


(4.5) 


(4.6) 


(4.7) 


The other equations needed for this system are derived by using the principle 
that the tangential derivatives of the shock conservation conditions across a wave 
are continuous and a theory developed by Lin and Rubinov in reference 21, which shows 
that if the curvature of a shock is known, then the spatial derivatives of the hydro- 
dynamical properties behind the shock can be found. 

Looking first at the weak wave f^, it has already been established that the 
properties u, v, p, T, and S are continuous across this wave for the case con- 
sidered. The tangential derivative along the wave x = f^(r) is 


-jj- ( )t. + -J- 

5x 4 dr 

r 

For each of the properties, then, the following equations hold: 


U 4 f 4 + U 4 = U 3 f 4 + U 3 
x r r x r r 


V. f A + V A = V 0 f A + V 0 

4 4 4 3 4 3 

x r r x r r 


P 4 f 4 + P 4 = P 3 f 4 + P 3 


(4.8) 

(4.9) 
(4.10) 


x r 


x r 


T. f . + T. = T, f . + T, 

4 4 4 3 4 3 

x r r x r r 


(4.11) 
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From the conditions at the slipstream 


Pa + P „ h ~ Prr + Pr- ^ 

4 4 x 5 5 x 

x r x r 


(4.1 2) 


and from the conditions of flow on either side of the slipstream 


u. h + u h - v. + u, h 2 - h =0 
4x 4 xx 4 4x 4x 

x x r r 


(4.1 3) 


u c h + u_h - v_ + u_ h 2 - v_. h =0 
5x 5 xx 5 5x 5x 
x x r r 


(4.14) 


The remaining equations in the system are derived from the following four shock 
equations: 


Continuity of mass 


p„T = p^T„ £> 

1 2 1 ^2 1 2 


r- Momentum 

(”s - S l) 

Energy 


*3 + V 5 - V 1 ‘ 0 

r 


(“s + “l)(“5 - 5 l) + ( ; 5 + ' 5 l) + (*5 + “l)(”5 - *l) 

2 T R(i 5 - *,) 


Y - 1 


= 0 


Combination of all shock equations 


2*yT 1 (l + f 2 + -| ) * (y + 1 )Q 1 Q 2 “ (Y ” 1 
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Taking the tangential derivative of each of these equations along the resultant 
shock f 3 yields the following four equations: 


P 1 T 5 2 1 - P 5 T 1 Q 5 - P 5 T 1 U 5 + P 5 T 1 V 5 f 3 f 3 
v x x x x r/ r 


+ P 1 T 5 2 1 ' P 5 T 1 2 5 - P 5 T 1 U 5 + P 5 T 1 V 5 f 3 
r r r r r 


+ PR T i f . + P i T * v i f ^ = I S 2 “ Pt “ Pi T 5 2 1 ] f 3 


513 ' ^1 A 5 V 1"3 

rr rr 


x r 


+ Vi e 5 - P 1 T 5 fi 1 - p ,Vl + P lVl f 3 + p lVl £ 

r r r r r 


rr 


(4.15) 


U 5 f 3 + V 5 ) f 3 + U 5 f 3 + U 5 f 3 + V 5 

{xr x r r r rr r 


- “t f 3 + V 1 f 3 + "l f 3 + "l f 3 + V 1 
x r xr rr rr r 


(4.16) 


V5 " 'Vi + V's "Vi + T^T^5 fs 

k x x x x ' x ' x r 


+ u„u„ - u u + v_v_ - V -V + 


5 5 1 1 5 5 1 1 y - 1 5 y - 

r r r r ' r ' 


T =0 (4.17) 


1 1 


2 Y 


'-T £ + T Wl + fM f 2* f f ] - f(Y * DIB, 0 5 

L x r r / V r / r rrj \ x 


+ Vs ' Q 1 V 5 f 3 - - D 

x x r / 


2Q 1 Q 1 "If 3 + |(y + D \ 

x r \ r 


“ V 1 2 5 f 3 2 5 V 1 f 3 + 2 1 U 5 2 1 V 5 f 3 " 2 1 V 5 f 3 

r r rr r r r rr/ 


(y - 1) 2 q7u - v f - v f 3 
\ r r r 


rr/ 


(4.18) 


32 



If it is assumed that all quantities in regions 1 and 3 are known, then the 

unknowns in the system are u,v,p,T,u,v,p, , T, 

44444444 

xxxxrrrr 

u 5 , v 5 , p 5 r T 5 , u 5 , v 5 , p 5 , T 5 r the curvature of the slipstream h xx , 

xxxxrrrr 

and the curvature of the resultant shock fo • *Ihis system of 1 8 linear equations 

rr 

is closed and has 18 unknowns, A matrix solution of this system yields the needed 
spatial derivatives for the system of equations (3.24) to (3.47). 


4.2 Spatial Derivatives in Regions 1 and 3 

For the spatial derivatives in regions 4 and 5 to be found, it is necessary that 

they be known in regions 1 and 3. By the use of characteristic equations and some 

difference methods, some of these have already been calculated in the MMOC program. 

In region 1, already available are u , v , and S . 

x x x 

Once it has been determined that an intersection occurs before the next C~ 
characteristic, then the characteristic subroutine is called as if the shock were not 
there. Data points 1 and 3 (see sketch F) are used to calculate the location of 



point 2 and all variables at point 2. A straight line interpolation is made along 
C“ between points 1 and 2 to get all values at Xe on that characteristic. Another 
interpolation is made along C* between points 2 and 3 to get values along this 
characteristic at Xe. Values at (Xe,Re) ahead of the first shock are then found by 
interpolating between Ra and Rb to get values at Re. Values of u r and v r 
are estimated in region 1 at (Xe,Re) by a simple difference equation between Ra and 
Rb. 


_ ua - ub 
u r “ Ra - Rb 


(4.19) 


- _ va - vb 
V r “ Ra - Rb 


(4.20) 
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With these values of u r and v r and the values of u x , v x , and S x already 
provided by the characteristic program, other spatial derivatives in region 1 are 
calculated in the following manner. From the x-momentum equation. 


p i 

x 



(J-Pi'Vi 




From the entropy equation, 


T 


1 


x 


( y - 

Y 




From the r-momentum equation. 


(4.21) 


( 4 . 22 ) 


P 

1 

r 



p u v 
1 1 1 


X 


p V V 
1 1 1 


r 


and from the continuity equation, 


i 'rrfiVi + Vi v i 

r P 1 v 1 \ x r 

... ... ViM 

+ V 1 T 1P1 - + — ; — / 

r x 


T p w 

1 


+ Vl P 1 


(4.23) 


(4.24) 


Equations (4.21) to (4.24) along with the derivatives already provided define all 
spatial derivatives needed in region 1 to determine the solution in regions 4 and 5. 


The spatial derivatives u, u, v, p, p, T, and T are also needed 

x r x x r x r 

in region 3 to provide information for the solution in regions 4 and 5. The values 
for fi , u , v^, v r , and are provided by the MMOC program using the char- 

acteristic equations behind the second shock. (See sketch G. ) 
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Again, the flow equations are used to provide the additional spatial derivatives 
needed in region 3 


x- Mo men turn 


x 


-P 3 U 3 U 3 " P 3 V 3 U 3 


(4.25) 


r-Momentum 


r RT 


" gp 3 " P 3 U 3 V 3 ' P 3 V 3 V 3 


(4.26) 


Entropy 


(Y + DT 3 /_ P 3 ? 


X 


Y i 3 + 


x p. 


(4.27) 


Continuity 


[T^p^U-x + T^p-V^ + 


T 3 P 3 W 3 


3 p v 3^3 3 ' 3^3 3 r 

r *3 3\ x r 


T 3 P 3 V 3 


+ U 3 T 3 P 3 + V 3 T 3 P 3 “ U 3 P 3 T 3 + r 


(4.28) 


Equations (4.25) to (4.28) and the derivatives provided by the characteristics 
in the MMOC program are used to determine the required spatial derivatives in 
region 3. There is now sufficient information available for the linear system in 
regions 4 and 5 to define all 16 spatial derivatives needed there. 


5. PROCEDURE FOR COMBINING COALESCENCE WITH MMOC PROGRAM 
5.1 Establishing Point of Coalescence 

In the MMOC program, calculations are performed from the bow shock along the C~ 
characteristic until either the last data point on the previous characteristic is 
reached or until the initial data line is intersected. (See sketch H.) The embedded 
shocks are either included as input, or else they occur during computation when char- 
acteristics of the same family cross each other. The program keeps track of the 
location of each embedded shock, and when each is reached, a separate routine is used 
to calculate the solution. The embedded shock routine uses the shock equations and 
the compatibility equation along the C + characteristic from the previous C” 
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Bow shock 


Embedded shocks 



C~(I) 


C (I-l) 


Sketch H 

characteristic to determine the solution just behind the shock. (See sketch I.) 
These equations determine the shock position, slope, curvature, and all variables 
behind the shock. 



The program monitors the location of all shocks, and when they are close to each 
other, a test is made to see if they cross each other before the next characteristic 
line. This test involves a straight line extrapolation of both shocks to estimate 
their intersection point Xe and an estimate of the point XI at which the first 
shock crosses the next C“ characteristic as shown in sketch J. If Xe is less 
than XI, control is passed to the coalescence routine with all data on the 1-1 th 
characteristic known and the data at points 1, 2, and 3 held. 
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c (I) 


Sketch J 

The following steps are then taken to determine the intersection point (Xe f Re) 
(see sketch K) : 


I) 


( 1 - 1 ) 


Sketch K 

(1) Interpolate between data points 1 , 2, and 3 to get values for all variables 
just ahead of the estimated intersection point in region 1. 

(2) Use the results of step 1 as data just ahead of the shock. Call the subrou- 

tine of the MMOC program which finds embedded shock solutions. This subroutine gives 
a local value of at Xe and the values of all the variables just behind the 

first shock. 

(3) Take the values of all variables behind shock 1 as the values in front of 
shock 2. Call the embedded shock subroutine to get the shock solution. This gives a 
value of P 2 and the values of all variables in region 3 at the intersection point. 
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(4) Adjust fa by taking p 1fnew = (fa + pi)/2 and p 2/ new = (p2 + Psi)/ 2 - 
Calculate a new value of (Xe,Re) using the adjusted values of p-j and P 2 . 

(5) Repeat steps 1 to 4 until the adjusted values of (3 change less than 

some error criterion e. Generally it is found that only three or four iterations 
are necessary with the error criterion on the order of 1CT 10 . 

When the above iteration has converged, the location of the coalescence point 
(Xe,Re) has been determined. Also known are the values of the variables in regions 1 
and 3, the local shock angles (3 1 and p 2 , and the local curvatures f 1 and f 2 at 
the intersection point. 


5.2 Local Solution 

Since the point of intersection and the data immediately ahead of this point are 
known, the spatial derivatives in regions 1 and 3 can now be found by the method 
described in section 4.2. The axisymmetric solution described in section 3.1 is now 
performed. This solution provides values for the variables u, v, p, T, and S 
in regions 4 and 5. It also provides the local slopes of f^, f^, and the slope of 
the slipstream h. The axisymmetric solution is independent of any of the tilde 
variables and w, and is found without iteration once the position of intersection is 
solved. The linear system described in section 4.1 is solved for the spatial deriva- 
tives of the variables in regions 4 and 5. There is now sufficient information to 
perform the asymmetric solution as described in section 3.2. This solution provides 
the values for u, v, p, T, and S in regions 4 and 5 for h x , h, f 3 , f 4 , 
and f 3 r » The local solution is now coirplete. 


5.3 Incorporation of Solution into Characteris tic Network 

The resultant shock f 3 is extrapolated linearly upward to obtain a first 
approximation of its intersection point with the Ith C“ characteristic, as shown in 
sketch L. Data at points 1 and 2 are known from the calculation described in 
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section 5,1 • A linear interpolation is used to find the value of all variables ahead 
of Xp, the intersection point. From the shock equations, the data ahead of the 
shock at Xp, and the C+ characteristic running from behind the shock at Xp to 
the top of the slipstream, a solution is found for £3 at Xp. new is calcu- 
lated using (($3 + (33)72, and the process is repeated until does not change 

between iterations. ' 

The point Xp and the data immediately behind the shock at Xp have now been 
determined. To retain the effects of the data below the slipstream, a new C“ char- 
acteristic (see sketch M) is calculated from the intersection point (Xe,Re) using the 
data from region 4. Ohe characteristic subroutine is called with data at points 1 



Sketch M 

and 2 to calculate point 3 , at points 1 * and 3 to calculate point 4 , and so on, until 
the end of the data on the 1-1 th characteristic is reached. 

At this point, it is recognized that the difference in entropy across the slip- 
stream is extremely small. It is convenient, therefore, to disregard the actual 
location of the slipstream and to allow the tangential discontinuity to become 
smeared in the region between C" ' and fg. The remainder of the Ith characteristic 
is now calculated using data from the intermediate characteristic C” . The effects 
of the data above the slipstream are felt in the location of Xp and in the data 
immediately behind Xp. (See sketch N. ) The effects of the data below the slip- 
stream appear 



Sketch N 
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in the value of the variables on the C~ characteris tic. To complete the Ith char- 
acteristic, the characteristic subroutine is called using points 1 and 2 to calculate 
point 3, and so on, until the characteristic is completed. Control of the program 
now reverts from the shock coalescence routine to the main program. If, instead of 
two embedded shocks intersecting, the bow shock intersects an embedded shock, the 
same procedure is followed with simplifications made because of the constant values 
ahead of the bow shock . 


6. RESULTS AND DISCUSSION 

The MMOC is a propagation method rather than a prediction method. Thus, flow- 
field input data on a cylindrical surface at approximately one body length from the 
axis must be supplied either from experiment or by use of computational methods. To 
employ the MMOC in an axisymmetric case, the initial data must include x, r, 0, 

M, and S/R at points on the cylinder in the plane (\> = 0. In addition, the 

strength of the bow shock and the location and strength of any embedded shocks must 
be given. For asymmetric cases, u, v, w, and S/R are also required. 

At the outset of the current work, it was felt that the flow-field computational 
method described in reference 25 represented the most accurate means available for 
predicting the data required as input to the MMOC. There also exist two alternate 
methods which can be used to calculate the field in the vicinity of the body, and 
both of these are also able to handle shock coalescence if they are allowed to calcu- 
late to the far field. The first, referred to as the modified uniform atmosphere 
method (MUAM) , is described in reference 26. In this method, based on Whitham the- 
ory, a uniform atmosphere is assumed, and a geometric mean correction for ambient 
pressure ^p a p g (where p^ is the pressure at the calculation altitude, and p^ is 
the ground pressure) is employed. In the second method, denoted as ARAP, geometric 
acoustics theory is used, and variable density in the atmosphere is accounted for. 
Each of these methods requires as input a description of the body and the flight 
conditions (ref. 27) . 

Since the axisymmetric shock coalescence is a new feature of the MMOC, it was 
initially decided to exercise this version of the program. Input data were provided 
for the MMOC at d/£ =0.2 by the method of reference 25 on a body of revolution to 
compare the extrapolated results of the MMOC with results predicted by the MUAM and 
ARAP programs. These axisymmetric results are discussed later in this section. 

The same body of revolution was placed at a small, positive angle of attack in 
order to obtain nonaxisymmetric data, and the method of reference 25 again was 
applied to calculate the required quantities at d/J£ = 0.2. At this point, two dif- 
ficulties arose which ultimately led to a decision to abandon the method of refer- 
ence 25 in favor of a less expensive, and presumably somewhat less accurate, determi- 
nation of the required data. First, the scheme of reference 25 is used to calculate 
quantities on a cylinder coaxial with the body and not on constant r, as required in 
the current work. Second, the presence of cross-derivative terms proportional to 
1/r throughout the asymmetric version of MMOC equations gives rise to numerical 
difficulties unless the initial surface is at least one body length away from the 
axis, a distance which requires excessive computation time when the method of refer- 
ence 25 is used for flow Mach numbers of 3 or 4. Since the aim of the current numer- 
ical computations was only to demonstrate the results of the new analysis, it was 
felt that continued use of this method for obtaining input data was not justified. 
Rather, it was decided to make use of the MUAM to determine appropriate data as input 
to the MMOC. The MUAM was modified to calculate cross-flow components, and the 
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derivatives with respect to c[> needed as input to the MMOC were calculated from 
their analytical expressions in the MUAM program. Accordingly, figure 1 includes the 
body used and the results obtained with the method of reference 25 for the axisymmet- 
ric case. All other results were obtained using a different body, and the MUAM was 
used to provide initial data for the MMOC at d/SL = 1.0. 

Axisymmetric data for the body of revolution shown in figure 1 were obtained 
using the method of reference 25 on a cylindrical surface at d /£ = 0.2* The pres- 
sure signature in the plane <J> = 0 at this radial location is shown. The MMOC pro- 
gram was then used to extrapolate this field through 900 body lengths to the ground. 
The resulting pressure signature is illustrated at the bottom of figure 1. It is 
emphasized that the evolution of the signature in figure 1 involves the coalescence 
of three shocks. Thus, even though it represents an axisymmetric case, this result 
cannot be obtained by using the MMOC without the analysis reported in the current 
paper. No comparison of this result with predictions of the MUAM or the ARAP pro- 
grams was made, since it was found that these methods had difficulty converging to a 
solution at 0.2 body lengths from the body. 

Correlations of MUAM and ARAP results with experimental data are shown in fig- 
ures 2 to 4. Figure 2 shows a comparison made in reference 26 for a body of revolu- 
tion at an angle of attack of 0° and a Mach number of 2.96. Though the bow shock 
amplitude is not predicted accurately, the expansion and tail pressures are well 
defined by the MUAM. Figure 3 shows a comparison of MUAM results with experimental 
data from reference 28. These data include angles of attack of 3° and -3° at Mach 
numbers of 1.41 and 2.01. The data points were read with pressure orifices on a 
reflection plate, and thus, a reflection factor of 2 has been used in the MUAM 
results. The experimental data show some scatter in these figures, but the pressure 
levels and trends are still rather well predicted by the MUAM at these distances for 
all angles of attack shown. No pressure data were available to verify pressures from 
the MUAM at an azimuthal angle other than 4> = 0. 

A correlation of flight-test data with ARAP predictions was made in refer- 
ence 7. Figure 4 is a reproduction of one such correlation made for Mach numbers 
of 1.35 to 3.00 over an altitude ranging from 35 000 feet to 70 000 feet. This 
figure shows that both the amplitude and the duration of the flight-test signatures 
were well predicted by the ARAP program. 

As originally written, the ARAP program was a sonic boom propagation method 
which required the Whitham F function as input. It was later modified to accept an 
area distribution which was transformed to the Whitham F function. This modifica- 
tion is based on having a smooth body with no discontinuities in slope. Though 
restricted to bodies of revolution, the MUAM program is not limited to smooth bodies. 
A second difference between the MUAM and the ARAP programs is that signals begin from 
the axis in the ARAP program and from the body surface in the MUAM program. Thus, 
when discontinuities are present, near-field shocks are more accurately located by 
the MUAM. 

With the experimental data shown then, it is felt that confidence can be placed 
in the pressure signature predictions by the MUAM in the near field for a body of 
revolution and in the ground level predictions made by the ARAP. Because there are 
no experimental data by which we can verify the results of the MMOC, comparisons are 
made of the MMOC results and the MUAM and ARAP results. Results from the MUAM at one 
body length were used to provide initial data for the MMOC. Comparisons of all three 
methods at several intermediate steps of propagation are shown. 
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The cone of revolution shown in figure 5 was used for the rest of the results. 

It has slope discontinuities at two locations so that embedded shocks are formed 
immediately, and thus, the shock coalescence portion of the program will be exer- 
cised. The pressure fields at one body length predicted by the MUAM and the ARAP 
for M = 3.00, a = 0° are shown at the upper left of figure 6. 

The two signatures are quite dissimilar. Recall that it is the MUAM signature 

which is felt to be valid in the near field and which is, therefore, used to initial- 
ize the MMOC. The remaining portions of figure 6 show comparisons of signatures 
predicted by all three methods during propagation. Note that at this Mach number and 
angle of attack, the MMOC and MUAM signatures remain very near each other through 
100 body lengths. At ground level, the bow shock predictions of the MMOC and the 

ARAP are very near each other, although the MMOC signature is predicted to be signif- 

icantly shorter. The positive and negative areas for the MMOC signature do not 
appear balanced in the ground signature in figure 6. However, the wave has been 
truncated, and considerable negative area has not been included. 

Asymmetric lifting forces are added through an angle of attack of 3° in fig- 
ure 7. Again, note that the MMOC and MUAM predictions remain very near each other in 
the near field (out to two body lengths), with the ARAP signature again being totally 
different. In the far-field and ground signatures, however, the shock amplitudes of 
the MMOC and ARAP results are very nearly equal, and the overall lengths of the sig- 
natures are not very different. 

It is noted here that comparisons are made of the shape and duration of the 
signature and not of its location in the atmosphere as denoted by the starting loca- 
tion x q . The three different methods of propagation do follow different paths 
through the atmosphere as determined by the position of x Q . By observing the loca- 
tions shown in figure 7, note that shock coalescence of the embedded and bow shocks 
has occurred between 2 and 10 body lengths. The ground signatures shown in these 
results all include a reflection factor of 1.8. 

In figure 8, even more asymmetries have been added to the flow field by 
increasing the angle of attack to 7°. The same general trend is seen in the near- 
field signatures. The difference in coalescence rate of the three methods can be 
observed in this figure. At two body lengths, the embedded shocks have coalesced in 
the MMOC and MUAM signatures but not in the ARAP signature, and at two and one-half 
body lengths, all shocks have coalesced in the MUAM and MMOC signatures. The ground 
signature again shows good correlation between the bow shock level of all signatures 
but at this angle of attack, the MMOC signature is now the longest of the three. 

Figures 9 and 10 show a similar set of comparisons at a = 0° and Mach numbers 
of 3.50 and 4.00, respectively. At M = 3.50, the MMOC and ARAP ground signatures 
have the same bow shock level, but the MMOC signature is somewhat shorter. At 
M = 4.00, the MMOC and ARAP ground signatures are practically identical in both 
amplitude and length. In both figures, the MUAM ground signature predicts the lowest 
bow shock and the longest signature. 

The experimental comparisons shown earlier for the MUAM and ARAP predictions 
indicate that confidence can be placed in the MUAM near-field predictions and in the 
ground-level ARAP predictions. These predictions are made based on a description of 
the equivalent body of revolution and on the flight conditions. The results from the 
MMOC program, when initialized with near-field MUAM data, approach ARAP ground shock 
predictions for all cases shown and have nearly the same signature as the ARAP at 
M = 4.00. 
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The fact that the MMOC program is shown to begin with valid data in the near 
field and to extrapolate to valid data (by comparison with the ARAP results) in the 
far field makes this a powerful tool in experimental work where predictions are made 
by the extrapolation of near-field signals, Ihe MMOC program is also able to make 
sonic boom predictions for complex three-dimensional bodies when their flow fields 
are provided by other computational methods. 

A comparison of the ground signatures when the MMOC program is initialized at 
two and one-half body lengths and at one body length, is shown in figure 11. Ihe 
signature extrapolated from two and one-half body lengths is shorter and has a lower 
bow shock than the signature extrapolated from one body length. This difference is 
to be expected because at two and one-half body lengths in figure 8, the MUAM signa- 
ture is already shorter than the MMOC signature and has a somewhat lower bow shock. 

The effect of asymmetries in the case with an angle of attack of 7° is seen in 
figure 12. The identical flow-field conditions with no asymmetric effects result in 
a signature whose bow shock is approximately 7 percent higher and whose length is 
3 percent shorter than the signature in which the asymmetries are included. 

Because one of the primary benefits in a program such as the MMOC would be for 
the extrapolation of wind-tunnel data to ground-level predictions, the effect of 
slight deviations in the determined locations of embedded shocks is an important 
consideration. In figure 13, the locations of the embedded shocks were perturbed 
about 1 percent, and the perturbed data were extrapolated to the ground. The results 
for the perturbed signature indicated approximately 2-percent deviation in the shock 
strength and in the length of the signature. 

Because a method based on linear theory was used to provide the initial data in 
these results, it was possible to calculate the necessary cross derivatives analyti- 
cally within the MUAM program. In the course of this study, however, several meth- 
ods, including curve fits and finite difference approximations, were used to compute 
these derivatives numerically, and it was found that the MMOC program is quite sensi- 
tive to their values. Thus, when experimental data or computational results which do 
not include analytically determined derivatives of the flow quantities are used to 
provide input to the MMOC program, extreme care must be taken in determining values 
for the cross derivative data. 


7. CONCLUDING REMARKS 

A method for analyzing shock coalescence including asymmetric effects has been 
presented. This method is based on an extension of the axisymmetric (locally two- 
dimensional) solution. The asymmetric effects are introduced through an additional 
set of governing equations, which are derived by taking the second circumferential 
derivative of the standard shock equations in the plane of symmetry. This shock 
coalescence method is consistent with and has been combined with a nonlinear sonic 
boom propagation method developed at New York University (NYU). The original NYU 
program, based on the method of characteristics, is unable to solve shock coalescence 
and ceases extrapolation when shock intersection occurs. This is often the case when 
data are extrapolated from realistic aircraft configurations. The combined program 
referred to as the modified method of characteristics (MMOC) is capable of extrapo- 
lating pressure signatures which include embedded shocks from the near field (approx- 
imately one body length from the axis) of an aircraft to ground level. Initial flow- 
field data are required on a cylindrical surface at approximately one body length 
from the aircraft axis. Flow variables and their cross derivatives are then defined 
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in the plane of symmetry beneath the aircraft and used during the extrapolation pro- 
cess. The location and strengths of all shocks on the initial surface must also be 
included in the data. 

The MMOC program is potentially extremely useful for predicting sonic booms from 
configurations where near-field data are provided either experimentally or from flow- 
field computational methods. Unlike previous extrapolation methods, the MMOC program 
accounts for the variation in entropy and nonlinear effects present near the aircraft 
and is thus able to begin extrapolation nearer to the body than the previous require- 
ment of about three body lengths. A major benefit of this capability is that sonic 
boom models may now be built approximately 3 feet in length (a restriction based on 
tunnel size) rather than the previous 6 inches. The larger size will allow a more 
detailed, realistic configuration than was previously possible. 

Included in the paper are extrapolated results from a body of revolution at a 
small, positive angle of attack. This configuration was designed so that embedded 
shocks would be included in the initial data and so that the MMOC solution could be 
compared with existing sonic boom prediction methods. Comparisons of signatures 
predicted by the MMOC with those predicted by two methods based on modified linear 
theory show good agreement in regions of propagation where the linear methods have 
been experimentally verified. The coalescence method as described contains weak 
shock/expansion approximations to determine spatial derivatives behind the wave of 
the opposite family, which may be needed to obtain the complete field in the vicinity 
of a point where shocks intersect. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
September 26, 1983 
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Figure 10.- Comparison of pressure signatures at M = 4.00 and a = 0° . 

as shown in figure 5; X = 100 ft; alt = 50 000 ft. 
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